## This file saves out .tex files of regression tables created in regression_models.R ##
## Created by Meredith Dost and last run 8/20/2025 ##

# load in packages and set up formatting
library(modelsummary); options("modelsummary_format_numeric_latex"="plain")
library(kableExtra)
library(broom)

format_nobs <- function(x) {
  if (is.numeric(x)) {
    return(format(x, big.mark = ",", scientific = FALSE))
  }
  return(x)
  }

gm <- list(
  list("raw" = "nobs", "clean" = "N", "fmt" = format_nobs),
  list("raw" = "adj.r.squared", "clean" = "Adj R$^2$", "fmt" = 3),
  list("raw" = "r2.within", "clean" = "Within R$^2$", "fmt" = 3))

# set working directory
# setwd("")

# load in results
load("mods_output.RData")

## main text condensed table
modelsummary(list(turn.med.mod.100,turn.bothidx.mod.100,turn.med_full.mod.100,turn.bothidx_full.mod.100), cluster = "fips",
             stars=T, 
             gof_map = gm, statistic = 'conf.int',
             coef_map = c("negsum" = "Medicaid burden index","exp"="Expanded Medicaid",
                          "regidx"="Registration burden index","turnidx"="Turnout burden index"),
             output = "tab2_indices_negsum_maintext.tex")

## full results among border counties
modelsummary(list(turn.med.mod.100,turn.bothidx.mod.100,turn.med_full.mod.100,turn.bothidx_full.mod.100), cluster = "fips",
             stars=T, 
             gof_map = gm, statistic = 'conf.int',
             coef_map = c("negsum" = "Medicaid burden index","exp"= "Expanded Medicaid",
                          "regidx"="Registration burden index","turnidx"="Turnout burden index",
                          "higheligib"="High Med eligibility","exp:higheligib"="Expansion*High elig interaction",
                          "pctwhite"="% non-Hispanic white",
                          "pcths"="% high school or less","pct65plus"="% age 65+","l_medincome"="Log median income",
                          "l_cvap"="Log voting eligible pop","demshare"="Democratic vote share","swing"="Swing state",
                          "sen"="Senatorial election","gub"="Gubernatorial election","sen:gub"="Sen*Gov interaction"),
             output = "tabA1_indices_negsum_border_appendix.tex")

## full results among all counties
modelsummary(list(turn.med.mod,turn.bothidx.mod,turn.med_full.mod,turn.bothidx_full.mod), cluster = "fips",
             stars=T, 
             gof_map = gm, statistic = 'conf.int',
             coef_map = c("negsum" = "Medicaid burden index","exp"= "Expanded Medicaid",
                          "regidx"="Registration burden index","turnidx"="Turnout burden index",
                          "higheligib"="High Med eligibility","exp:higheligib"="Expansion*High elig interaction",
                          "pctwhite"="% non-Hispanic white",
                          "pcths"="% high school or less","pct65plus"="% age 65+","l_medincome"="Log median income",
                          "l_cvap"="Log voting eligible pop","demshare"="Democratic vote share","swing"="Swing state",
                          "sen"="Senatorial election","gub"="Gubernatorial election","sen:gub"="Sen*Gov interaction"),
             output = "tabA2_indices_negsum_allctys_appendix.tex")

## alternative Medicaid burden indices
modelsummary(list(turn.bothidx.mod.f,turn.bothidx_full.mod.f, turn.bothidx.mod.r3 ,turn.bothidx_full.mod.r3, 
                  turn.bothidx.mod.100.f,turn.bothidx_full.mod.100.f,
                  turn.bothidx.mod.100.r3, turn.bothidx_full.mod.100.r3), cluster = "fips",
             stars=T, 
             gof_map = gm, statistic = 'conf.int',
             coef_map = c("foxneg" = "Alternative Medicaid index",
                          "AB3" = "Reduced burden index","exp"= "Expanded Medicaid",
                          "regidx"="Registration burden index","turnidx"="Turnout burden index",
                          "higheligib"="High Med eligibility","exp:higheligib"="Expansion*High elig interaction",
                          "pctwhite"="% non-Hispanic white",
                          "pcths"="% high school or less","pct65plus"="% age 65+","l_medincome"="Log median income",
                          "l_cvap"="Log voting eligible pop","demshare"="Democratic vote share","swing"="Swing state",
                          "sen"="Senatorial election","gub"="Gubernatorial election","sen:gub"="Sen*Gov interaction"),
             output = "tabA3_indices_versions_FULL.tex")

## COVI & presidential election + midterm year models index
modelsummary(list(turn.bothidx_full.mod.cv,turn.bothidx_full.mod.p,turn.bothidx_full.mod.m,turn.bothidx_full.mod.100.cv,turn.bothidx_full.mod.100.p,turn.bothidx_full.mod.100.m), cluster = "fips",
             stars=T, 
             gof_map = gm, statistic = 'conf.int',
             coef_map = c("negsum" = "Medicaid burden index","exp"= "Expanded Medicaid",
                          "FinalCOVI"="COVI",
                          "regidx"="Registration burden index","turnidx"="Turnout burden index",
                          "higheligib"="High Med eligibility","exp:higheligib"="Expansion*High elig interaction",
                          "pctwhite"="% non-Hispanic white",
                          "pcths"="% high school or less","pct65plus"="% age 65+","l_medincome"="Log median income",
                          "l_cvap"="Log voting eligible pop","demshare"="Democratic vote share","swing"="Swing state",
                          "sen"="Senatorial election","gub"="Gubernatorial election","sen:gub"="Sen*Gov interaction"),
             output = "tabA4_covi_pres_mid_results.tex")

## combined election index
modelsummary(list(turn.bothidx.mod.e,turn.bothidx_full.mod.e,turn.bothidx.mod.100.e,turn.bothidx_full.mod.100.e), cluster = "fips",
             stars=T, 
             gof_map = gm, statistic = 'conf.int',
             coef_map = c("negsum" = "Medicaid burden index","exp"= "Expanded Medicaid",
                          "cov"="Combined election index",
                          "higheligib"="High Med eligibility","exp:higheligib"="Expansion*High elig interaction",
                          "pctwhite"="% non-Hispanic white",
                          "pcths"="% high school or less","pct65plus"="% age 65+","l_medincome"="Log median income",
                          "l_cvap"="Log voting eligible pop","demshare"="Democratic vote share","swing"="Swing state",
                          "sen"="Senatorial election","gub"="Gubernatorial election","sen:gub"="Sen*Gov interaction"),
             output = "tabA5_elec_indices_combined.tex")

## weighted results
modelsummary(list(turn.bothidx.mod.wt,turn.bothidx_full.mod.wt,turn.bothidx.mod.100.wt,turn.bothidx_full.mod.100.wt), cluster = "fips",
             stars=T, 
             gof_map = gm, statistic = 'conf.int',
             coef_map = c("negsum" = "Medicaid burden index","exp"= "Expanded Medicaid",
                          "regidx"="Registration burden index","turnidx"="Turnout burden index",
                          "higheligib"="High Med eligibility","exp:higheligib"="Expansion*High elig interaction",
                          "pctwhite"="% non-Hispanic white",
                          "pcths"="% high school or less","pct65plus"="% age 65+","l_medincome"="Log median income",
                          "l_cvap"="Log voting eligible pop","demshare"="Democratic vote share","swing"="Swing state",
                          "sen"="Senatorial election","gub"="Gubernatorial election","sen:gub"="Sen*Gov interaction"),
             output = "tabA8_weighted_results.tex")

## interaction model results
modelsummary(list(turn.bothidx.mod.mp, turn.bothidx_full.mod.mp, turn.bothidx.mod.100.mp, turn.bothidx_full.mod.100.mp), cluster = "fips",
             stars=T, 
             gof_map = gm, statistic = 'conf.int',
             coef_map = c("negsum" = "Medicaid burden index","exp"= "Expanded Medicaid",
                          "regidx"="Registration burden index","turnidx"="Turnout burden index",
                          "pcteligib"="% Med expansion eligible","negsum:pcteligib"="Medicaid burden*% eligible interaction",
                          "exp:higheligib"="Expansion*High elig interaction",
                          "pctwhite"="% non-Hispanic white",
                          "pcths"="% high school or less","pct65plus"="% age 65+","l_medincome"="Log median income",
                          "l_cvap"="Log voting eligible pop","demshare"="Democratic vote share","swing"="Swing state",
                          "sen"="Senatorial election","gub"="Gubernatorial election","sen:gub"="Sen*Gov interaction"),
             output = "tabA9_interaction_results.tex")

#cleaning up R environment
rm(list=ls()[!ls() %in% c("format_nobs","gm","turn.bothidx.mod.100","turn.bothidx_full.mod.100")])
################ load in bootstrapped models #########
load("bootstraps_mods_output.RData")
modelsummary(list(turn.bothidx.mod.100,turn.bothidx.mod.100,turn.bothidx_full.mod.100,turn.bothidx_full.mod.100), cluster = "fips",
             stars=T, 
             gof_map = gm, statistic = 'conf.int',
             coef_map = c("negsum" = "Medicaid burden index"),
             output = "tabA7_twfe_comparison.tex")

## note that you need to manually input into the above .tex file the coefficients and results that use the de Chaisemartin et al 2024 estimators.
round(mod.tot.nocovar$results$Effects,3)
Nm=2447+3352
t=mod.tot.nocovar$results$Effects[1]/mod.tot.nocovar$results$Effects[2]
dt(t,df=Nm-length(coef(turn.bothidx.mod.100))-2-2791)

round(mod.tot$results$Effects,3)
t=mod.tot$results$Effects[1]/mod.tot$results$Effects[2]
dt(t,df=Nm-length(coef(turn.bothidx_full.mod.100))-2-2791)

#cleaning up R environment
rm(list=ls()[!ls() %in% c("format_nobs")])
##### heterogeneous treatment effect diagnostics #######
load("heterogeneouschecks.RData")

modelsummary(check_slopes,
             coef_rename = c("treatment_resid" = "Residualized treatment",
                             "negsum" = "Medicaid burden index",
                             "treatment_resid:negsum" = "Medicaid burden index × residualized treatment",
                             "(Intercept)" = "Intercept"),
             estimate = "{estimate}",
             stars=T, 
            statistic = 'conf.int',
             gof_map = tribble(
               ~raw, ~clean, ~fmt,
               "nobs", "N", 0,
               "r.squared", "R²", 3),
            output = "tabA6_hetero_slopes.tex")


#cleaning up R environment
rm(list=ls()[!ls() %in% c("format_nobs")])
##### CCES analyses #######
load("cces_output.RData")
gm <- list(
  list("raw" = "nobs", "clean" = "N", "fmt" = format_nobs),
  list("raw" = "adj.r.squared", "clean" = "Adj R$^2$", "fmt" = 3),
  list("raw" = "logLik", "clean" = "Log Likelihood", "fmt" = 3),
  list("raw" = "bic", "clean" = "BIC", "fmt" = 3))

modelsummary(list(out_d2, out, out_not), exponentiate=F,
             stars=T, 
             gof_map = gm, statistic = 'conf.int',
             coef_map = c("negsum1612" = "Medicaid burden index",
                          "negsum1612:highmed" = "Med burden*High Med enrollment ", "highmed"="High Med enrollment",
                          "regidx"="Registration burden index","turnidx"="Turnout burden index",
                          "exp_altfirst"= "Expanded Medicaid in 2014","exp_altnext"= "Expanded Medicaid after 2014",
                          "higheligib" = "High Med eligibility",
                          "exp_altfirst:higheligib" = "Expans in 2014*High eligibility",
                          "exp_altnext:higheligib" = "Expans after 2014*High eligibility",
                          "blacknh"="Black, non-Hispanic","othernh"="Other race, non-Hispanic",
                          "hisp"="Hispanic","educ"="Education level","age"="Age","faminc"="Family income",
                          "female"="Female","med"="Medicaid recipient"),
             output = "tab3_cces_ecological_expo.tex")

## all adults, including 65+
modelsummary(list(out_d.all), exponentiate=F,
             stars=T, 
             gof_map = gm, statistic = 'conf.int',
             coef_map = c("negsum1612" = "Medicaid burden index",
                          "negsum1612:highmed" = "Med burden*High Med enrollment ", "highmed"="High Med enrollment",
                          "regidx"="Registration burden index","turnidx"="Turnout burden index",
                          "exp_altfirst"= "Expanded Medicaid in 2014","exp_altnext"= "Expanded Medicaid after 2014",
                          "higheligib" = "High Med eligibility",
                          "exp_altfirst:higheligib" = "Expans in 2014*High eligibility",
                          "exp_altnext:higheligib" = "Expans after 2014*High eligibility",
                          "blacknh"="Black, non-Hispanic","othernh"="Other race, non-Hispanic",
                          "hisp"="Hispanic","educ"="Education level","age"="Age","faminc"="Family income",
                          "female"="Female","med"="Medicaid recipient"),
             output = "tabA10_cces_ecological_totaladults_expo.tex")

